---
title: "regression"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
library(haven)
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)

theme_set(theme_sjplot())

dat <- read_dta("data/iraq2022_ver1.dta")
dat2019 <- read_dta("data/iraq2019.dta")

dat <- dat %>% 
  mutate(Security_rely = recode(security_rely,
                      "1" = "Govt",
                      "2" = "Military",
                      "3" = "Police",
                      "4" = "PMU",
                      "5" = "others"))
dat <- dat %>% 
  mutate(Security_rely = fct_relevel(Security_rely,
                                  "Govt", "Military", "Police", "PMU", "others"))

dat2019 <- dat2019 %>% 
  mutate(Security_rely = recode(security_rely,
                      "1" = "Govt",
                      "2" = "Military",
                      "3" = "Police",
                      "4" = "PMU",
                      "5" = "others"))

dat2019 <- dat2019 %>% 
  mutate(Security_rely = fct_relevel(Security_rely,
                                  "Govt", "Military", "Police", "PMU", "others"))
dat <- dat %>% 
  mutate(sect = as.numeric(D9))

dat <- dat %>% 
  mutate(Trust_police = recode(trust_police,
                      "1" = "not at all",
                      "2" = "not",
                      "3" = "nutral",
                      "4" = "trust",
                      "5" = "very much"))
dat <- dat %>% 
  mutate(Trust_police = fct_relevel(Trust_police,
                                  "not at all", "not", "nutral", "trust", "very much"))

```

# overview
```{r}
stargazer(as.data.frame(dat), type = "text")
stargazer(as.data.frame(dat2019), type = "text")
```

# trust Police
```{r}
barplot(table(dat$shia, dat$Trust_police),
        main = "Trust police",
        beside = TRUE)
legend("topleft", col = c("black", "gray"), pch = 15, legend = c("Shia","others"))   
```

# rely security on PMU
```{r}
par(mfrow = c(1, 2))

t1 <- barplot(table(dat2019$Shia, dat2019$Security_rely),
              main = "Rely in security in 2019 (N = 1150)",
              beside = TRUE)
legend("topright", col = c("black", "gray"), pch = 15, legend = c("Shia","others"))   

t2 <- barplot(table(dat$shia, dat$Security_rely),
        main = "Rely in security in 2022 (N = 1500)",
              beside = TRUE)
legend("topright", col = c("black", "gray"), pch = 15, legend = c("Shia","others"))   

```

# overview dependent variavles
```{r}
hist(dat$PMU_experiment)
```

# regression: shia/minorities
```{r}
reg01 <- lm(formula = PMU_experiment ~ list_PMU, data = dat)
summary(reg01)

reg02 <- lm(formula = PMU_experiment ~ list_PMU*shia, data = dat)
summary(reg02)

reg03 <- lm(formula = PMU_experiment ~ list_PMU*others, data = dat)
summary(reg03)

reg04 <- lm(formula = PMU_experiment ~ list_PMU*shia + 
              list_PMU*others, data = dat)
summary(reg04)

reg05 <- lm(formula = PMU_experiment ~ list_PMU*shia + 
              sex + age + education + income, data = dat)
summary(reg05)

reg06 <- lm(formula = PMU_experiment ~ list_PMU*others + 
              sex + age + education + income, data = dat)
summary(reg06)

reg07 <- lm(formula = PMU_experiment ~ list_PMU*sunni + 
              sex + age + education + income, data = dat)
summary(reg07)

reg08 <- lm(formula = PMU_experiment ~ list_PMU*shia + 
              list_PMU*sunni + list_PMU*others +
              sex + age + education + income, data = dat)
summary(reg08)



reg011 <- lm_robust(formula = PMU_experiment ~ list_PMU, 
                   cluster = C5, se_type = "stata",
                   data = dat)
summary(reg011) 

reg021 <- lm_robust(formula = PMU_experiment ~ list_PMU*shia, 
                   cluster = C5, se_type = "stata",
                   data = dat)
summary(reg021) 

reg031 <- lm_robust(formula = PMU_experiment ~ list_PMU*others, 
                   cluster = C5, se_type = "stata",
                   data = dat)
summary(reg031) 

reg041 <- lm_robust(formula = PMU_experiment ~ list_PMU*shia +
                      list_PMU*others, 
                   cluster = C5, se_type = "stata",
                   data = dat)
summary(reg041) 

reg051 <- lm_robust(formula = PMU_experiment ~ list_PMU*shia +
                      list_PMU*others +
                      sex + age + education + income, 
                   cluster = C5, se_type = "stata",
                   data = dat)
summary(reg051) 
```

# regression Iran/US
```{r}
reg20 <- lm(formula = PMU_experiment ~ list_PMU*Iran, data = dat)
summary(reg20)

reg21 <- lm(formula = PMU_experiment ~ list_PMU*US, data = dat)
summary(reg21)

reg22 <- lm(formula = PMU_experiment ~ list_PMU*Iran + 
              sex + age + education + income, data = dat)
summary(reg22)

reg23 <- lm(formula = PMU_experiment ~ list_PMU*US +
              sex + age + education + income, data = dat)
summary(reg23)
```

# regression Fatah/Sadr
```{r}
reg30 <- lm(formula = PMU_experiment ~ list_PMU*Fatah, data = dat)
summary(reg30)

reg31 <- lm(formula = PMU_experiment ~ list_PMU*Sadr, data = dat)
summary(reg31)

reg32 <- lm(formula = PMU_experiment ~ list_PMU*Fatah + 
              list_PMU*Sadr, data = dat)
summary(reg32)

reg33 <- lm(formula = PMU_experiment ~ list_PMU*Fatah + 
              sex + age + education + income, data = dat)
summary(reg33)
```

# regression trust
```{r}
reg40 <- lm(formula = PMU_experiment ~ list_PMU*trust_military, data = dat)
summary(reg40)

reg41 <- lm(formula = PMU_experiment ~ list_PMU*trust_police, data = dat)
summary(reg41)

reg42 <- lm(formula = PMU_experiment ~ list_PMU*trust_religious_leader, data = dat)
summary(reg42)

reg43 <- lm(formula = PMU_experiment ~ list_PMU*trust_military + 
              list_PMU*trust_police + list_PMU*trust_religious_leader, data = dat)
summary(reg43)

reg44 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income, data = dat)
summary(reg44)

# shia
reg45 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, sect == "2"))
summary(reg45)

# sunni
reg46 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, sect == "1"))
summary(reg46)

# kurd
reg47 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, sect == "4"))
summary(reg47)

```

# regression rely
```{r}
reg50 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_robbery, data = dat)
summary(reg50)

reg51 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_job, data = dat)
summary(reg51)

reg52 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security, data = dat)
summary(reg52)

reg53 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, data = dat)
summary(reg53)

# shia
reg54 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, sect == "2"))
summary(reg54)

# sunni
reg55 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, sect == "1"))
summary(reg55)

# Kurd
reg56 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, sect == "4"))
summary(reg56)
```

# regression party
```{r}
reg60 <- lm(formula = PMU_experiment ~ list_PMU*support_shiaparties +
              sex + age + education + income, data = dat)
summary(reg60)

reg61 <- lm(formula = PMU_experiment ~ list_PMU*support_sunniparties +
              sex + age + education + income, data = dat)
summary(reg61)

reg62 <- lm(formula = PMU_experiment ~ list_PMU*vote_shiaparties +
              sex + age + education + income, data = dat)
summary(reg62)

reg63 <- lm(formula = PMU_experiment ~ list_PMU*vote_sunniparties +
              sex + age + education + income, data = dat)
summary(reg63)

reg64 <- lm(formula = PMU_experiment ~ list_PMU*no_party_support +
              sex + age + education + income, data = dat)
summary(reg64)
```


# regression region
```{r}
reg70 <- lm(formula = PMU_experiment ~ list_PMU*North +
              sex + age + education + income, data = dat)
summary(reg70)

reg71 <- lm(formula = PMU_experiment ~ list_PMU*Central +
              sex + age + education + income, data = dat)
summary(reg71)

reg72 <- lm(formula = PMU_experiment ~ list_PMU*West +
              sex + age + education + income, data = dat)
summary(reg72)

reg73 <- lm(formula = PMU_experiment ~ list_PMU*South +
              sex + age + education + income, data = dat)
summary(reg73)

reg74 <- lm(formula = PMU_experiment ~ list_PMU*KRG +
              sex + age + education + income, data = dat)
summary(reg74)

reg700 <- lm(formula = PMU_experiment ~ list_PMU*North +
               list_PMU*Central +
               list_PMU*West +
               list_PMU*South +
              sex + age + education + income, data = dat)
summary(reg700)

# rely on security
# north
reg75 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, region == "2"))
summary(reg75)

# central
reg76 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, region == "3"))
summary(reg76)

# west
reg77 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, region == "4"))
summary(reg77)

# south
reg78 <- lm(formula = PMU_experiment ~ list_PMU*rely_PMU_security + 
              sex + age + education + income, 
            subset(dat, region == "5"))
summary(reg78)


# trust police
# north
reg80 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, region == "2"))
summary(reg80)

# central
reg81 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, region == "3"))
summary(reg81)

# west
reg82 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, region == "4"))
summary(reg82)

# south
reg83 <- lm(formula = PMU_experiment ~ list_PMU*trust_police + 
              sex + age + education + income,
            subset(dat, region == "5"))
summary(reg83)
```


# regression table
```{r}
reg_table <- mtable("Treatment" = reg01, "Sect" = reg08,
                    "Support Shia parties" = reg60, "Fatah" = reg33, 
                    "Intervension" = reg22, 
               summary.stats=c("adj. R-squared","F","p", "N"))
print(reg_table, type = "html")
write_html(reg_table, file = "result.html")
write.mtable(reg_table, file = "result.csv", colsep = ",")

reg_region <- mtable("North" = reg70, "Central" = reg71, "West" = reg72, 
                    "South" = reg73, "LRG" = reg74,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(reg_region, type = "html")
write_html(reg_region, file = "result_resion.html")
write.mtable(reg_region, file = "result_region.csv", colsep = ",")
```

# coef_plot
```{r}
multiplot(reg01, reg05, reg60, reg33, reg22, reg44, reg53,
          title = "OLS results",
          names = c("M1:treatment", "M2:sect", "M3:party", "M4:Fatah", "M5:external", "M6:trust", "M7:rely"),
          sort = c("magnitude"),
          intercept = FALSE, plot.linetypes = FALSE, legend.reverse = FALSE,
          plot.shapes = FALSE,
          numberAngle = 0, lwdOuter = 1, 
          coefficients = c("list_PMU", 
                           "list_PMU:shia",
                           "list_PMU:support_shiaparties",
                           "list_PMU:Fatah",
                           "list_PMU:Iran",
                           "list_PMU:trust_police",
                           "list_PMU:rely_PMU_security"),
          newNames = c(list_PMU = "Treatment group", 
                       "list_PMU:shia" = "Shia",
                       "list_PMU:support_shiaparties" = "Support Shia parties",
                       "list_PMU:Fatah" = "Support Fatah",
                       "list_PMU:Iran" = "Iranian intervention",
                       "list_PMU:trust_police" = "Trust police",
                       "list_PMU:rely_PMU_security" = "Rely PMU"),
          single = FALSE,
          ncol = 4) +
  theme_bw() +
  theme(legend.position = "none")+
  ggtitle("")
```


# marginal effect: sect
```{r}
int_reg01 <- interplot(m = reg05, var1 = "list_PMU", var2 = "shia", point = TRUE) +
  labs(x = "Shia ***",
       y = "Effect of treatment") +
  ggtitle("Model 2 (Shia)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg01)

p1 <- plot_model(reg05, type = "pred", terms = c("shia", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(sect)")

int_reg011 <- interplot(m = reg06, var1 = "list_PMU", var2 = "others", point = TRUE) +
  labs(x = "Minority †",
       y = "Effect of treatment") +
  ggtitle("Model 2 (Minority)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg011)

p11 <- plot_model(reg06, type = "pred", terms = c("others", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(sect)")

int_reg012 <- interplot(m = reg07, var1 = "list_PMU", var2 = "sunni", point = TRUE) +
  labs(x = "Sunni",
       y = "Effect of treatment") +
  ggtitle("Model 2 (Sunni)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg012)

p12 <- plot_model(reg07, type = "pred", terms = c("sunni", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(sect)")
```

# marginal effect: external
```{r}
int_reg02 <- interplot(m = reg22, var1 = "list_PMU", var2 = "Iran", point = TRUE) +
  labs(x = "Iranian intervention",
       y = "Effect of treatment") +
  ggtitle("Model 5") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg02)

p2 <- plot_model(reg22, type = "pred", terms = c("Iran", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(Iranian intervention)")
```

# marginal effect: Fatah
```{r}
int_reg03 <- interplot(m = reg33, var1 = "list_PMU", var2 = "Fatah", point = TRUE) +
  labs(x = "Support Fatah*",
       y = "Effect of treatment") +
  ggtitle("Model 4") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg03)

p3 <- plot_model(reg33, type = "pred", terms = c("Fatah", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(support Fatah)")

reg330 <- lm(formula = PMU_experiment ~ Fatah*list_PMU + 
              sex + age + education + income, data = dat)
summary(reg330)

plot_model(reg330, type = "int",
           axis.labels = "Treatment", title = "Subs. eff.(support Fatah)")

plot_model(reg33, type = "pred", terms = c("Fatah", "list_PMU", "income"),
           axis.labels = "Treatment", title = "Subs. eff.(support Fatah)")
```

# marginal effect: trust
```{r}
int_reg04 <- interplot(m = reg44, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, overall***)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg04)

p4 <- plot_model(reg44, type = "pred", terms = c("trust_police", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(trust police)")


# shia
int_reg041 <- interplot(m = reg45, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, Shia*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg041)

# sunni
int_reg042 <- interplot(m = reg46, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, Sunni*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg042)

# kurd
int_reg043 <- interplot(m = reg47, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, Kurd)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg043)


# north
int_reg044 <- interplot(m = reg80, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, North)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg044)

# central
int_reg045 <- interplot(m = reg81, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, Central***)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg045)

# west
int_reg046 <- interplot(m = reg82, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, West*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg046)

# south
int_reg047 <- interplot(m = reg83, var1 = "list_PMU", var2 = "trust_police", point = TRUE) +
  labs(x = "Trust police",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(trust police, South)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg047)
```

# marginal effect: rely
```{r}
int_reg05 <- interplot(m = reg53, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, overall*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg05)

p5 <- plot_model(reg53, type = "pred", terms = c("rely_PMU_security", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(rely PMU)")

# shia
int_reg051 <- interplot(m = reg54, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, Shia*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg051)

# sunni
int_reg052 <- interplot(m = reg55, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, Sunni*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg052)

# Kurd
int_reg053 <- interplot(m = reg56, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, Kurd)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg053)


# north
int_reg054 <- interplot(m = reg75, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, North***)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg054)

# central
int_reg055 <- interplot(m = reg76, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, Central*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg055)

# west
int_reg056 <- interplot(m = reg77, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, West)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg056)

# sounth
int_reg057 <- interplot(m = reg78, var1 = "list_PMU", var2 = "rely_PMU_security", point = TRUE) +
  labs(x = "Rely on PMU in security",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(rely PMU, South*)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg057)
```

# marginal effect: party
```{r}
int_reg06 <- interplot(m = reg60, var1 = "list_PMU", var2 = "support_shiaparties", point = TRUE) +
  labs(x = "Support Shia parties***",
       y = "Effect of treatment") +
  ggtitle("Model 3") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg06)

p6 <- plot_model(reg60, type = "pred", terms = c("support_shiaparties", "list_PMU"), 
           axis.labels = "Treatment", title = "Subs. eff.(support Shia parties)")
```

# marginal effect: region
```{r}
int_reg07 <- interplot(m = reg70, var1 = "list_PMU", var2 = "North", point = TRUE) +
  labs(x = "North*",
       y = "Effect of treatment") +
  ggtitle("Model 6") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg07)

int_reg08 <- interplot(m = reg71, var1 = "list_PMU", var2 = "Central", point = TRUE) +
  labs(x = "Central**",
       y = "Effect of treatment") +
  ggtitle("Model 7") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg08)

int_reg09 <- interplot(m = reg72, var1 = "list_PMU", var2 = "West", point = TRUE) +
  labs(x = "West**",
       y = "Effect of treatment") +
  ggtitle("Model 8") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg09)

int_reg10 <- interplot(m = reg73, var1 = "list_PMU", var2 = "South", point = TRUE) +
  labs(x = "South",
       y = "Effect of treatment") +
  ggtitle("Model 9") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg10)

int_reg11 <- interplot(m = reg74, var1 = "list_PMU", var2 = "KRG", point = TRUE) +
  labs(x = "KRG",
       y = "Effect of treatment") +
  ggtitle("Model 10") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg11)
```

# plot marginal effect
```{r}
ggarrange(int_reg01, int_reg012, int_reg011,
          int_reg06, int_reg03, int_reg02,
          nrow = 2, ncol = 3, align = "hv")
```

# plot marginal effect, region
```{r}
ggarrange(int_reg07, int_reg08, int_reg09,
          int_reg10, int_reg11,
          nrow = 2, ncol = 3, align = "hv")
```

# plot marginal effect, trust police, sect
```{r}
ggarrange(int_reg04, int_reg041, 
          int_reg042, int_reg043,
          nrow = 2, ncol = 2, align = "hv")
```

# plot marginal effect, trust police, region
```{r}
ggarrange(int_reg044, int_reg045, 
          int_reg046, int_reg047,
          nrow = 2, ncol = 2, align = "hv")
```

# plot marginal effect, rely on, sect
```{r}
ggarrange(int_reg05, int_reg051, 
          int_reg052, int_reg053,
          nrow = 2, ncol = 2, align = "hv")
```

# plot marginal effect, rely on, region
```{r}
ggarrange(int_reg054, int_reg055, 
          int_reg056, int_reg057,
          nrow = 2, ncol = 2, align = "hv")
```

# plot substantial effect
```{r}
ggarrange(p1, p6, p3,
          p2, p4, p5,
          nrow = 2, ncol = 3, align = "hv")
```

# appendix
```{r}
ap_dat <- dat %>% 
  select(list_PMU, shia, sunni, others, sex, age, education, income,
         Fatah, Iran, North, Central, West, South, KRG)
namelist <- names(ap_dat)
ap_dat[, namelist] <- lapply(ap_dat[, namelist], as.numeric)
stargazer(as.data.frame(ap_dat), type = "text")
```


# H1 robust: ordinal regression
```{r}
library(MASS)
library(effects)

reg61 <- polr(as.factor(PMU_experiment) ~ list_PMU, data = dat)
summary(reg61)

reg62 <- polr(as.factor(PMU_experiment) ~ list_PMU*others +
            sex + age + education + income, data = dat)
summary(reg62)
(ctable <- coef(summary(reg62)))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE)*2
(ctable <- cbind(ctable, "p value" = p))

library(ordinal)
reg71 <- clm(as.factor(PMU_experiment) ~ list_PMU*Fatah, data = dat)
summary(reg71)

reg72 <- clm(as.factor(PMU_experiment) ~ list_PMU*trust_police +
            sex + age + education + income, data = dat)
summary(reg72)
```


